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Abstract. There exists, in general, a convex set of quantum state estimators 
that maximize the likelihood for informationally incomplete data. We propose 
an estimation scheme, catered to measurement data of this kind, to search for 
the exact maximum-likelihood-maximum-entropy estimator using semidcfinitc 
programming and a standard multi-dimensional function optimization routine. 
This scheme can be used to infer the expectation values of a set of entanglement 
witnesses that can be used to verify the entanglement of the unknown quantum 
state for composite systems. Next, we establish an alternative numerical scheme 
that is more computationally robust for the sole purpose of maximizing the 
likelihood and entropy. 
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1. Introduction 

Quantum state preparation is the first important step for any protocol that makes use 
of quantum resources. Examples of such protocols are quantum state teleportation 
and quantum key distribution which require entangled quantum states. In order to 
verify the integrity of the true quantum state ptruo prepared by the source, one carries 
out quantum state tomography to characterize it. Measurements are performed on a 
collection of identical copies of quantum systems (electrons, photons, etc.) that are 
emitted from the source. Then, the quantum state of the source is inferred from the 
measurement data obtained from this collection. The measurements are generically 
described by a set of positive operators {lb,} that compose a probability operator 
measurement (POM). Such a procedure of state inference is known as quantum state 
estimation. 

When the measurement outcomes form an informationally complete set, they fully 
characterize the source and the measurement data obtained will contain maximal 
information about its state. To infer the unknown state from the data, one can 
search for a state estimator that maximizes the likelihood functional that yields 
the probability of obtaining a particular sequence of measurement detections given 
a quantum state - the maximum-likelihood (ML) estimator [TJ El HJ ■ Yet, in 
tomography experiments performed on complex quantum systems with many degrees 
of freedom, it is not possible to implement such an informationally complete set 
of measurement outcomes. Therefore, some information about the source will be 
missing and its quantum state cannot be unambiguously determined. For instance, if 
a source produces a mode of light that is described by an infinite-dimensional statistical 
operator ptrue, then no matter how ingeniously a measurement scheme is designed to 
probe incoming photons prepared by this source, an infinite amount of information 
about the mode of light will always remain unknown. The ML estimator obtained from 
these informationally incomplete data is no longer unique and there will in general 
be infinitely many other ML estimators that are consistent with the data. These 
estimators form a convex set under the likelihood plateau. 

In order to choose an estimator from the convex set for statistical prediction, we 
can consider the maximum-entropy principle advocated by E. T. Jaynes [3 |6]. In 
doing so, one obtains a unique estimator that maximizes both the likelihood and 
the von Neumann entropy functional. Statistically, this estimator is least-biased 
for the informationally incomplete data. In E], we developed and applied 
an algorithm, based on the steepest-ascent method, to approximately look for the 
maximum-likelihood-maximum-entropy (MLME) estimator. This algorithm involves 
a parameter that needs to be chosen just above a minimum threshold to obtain an 
estimator that is as close to the actual MLME estimator as possible. In general, this 
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threshold depends on the true state ptmc- Therefore, one needs to run the algorithm 
a few times to estimate this threshold. 

This paper is organized as follows. After a brief review of this steepest-ascent 
algorithm in Section [3J we introduce a numerical scheme that is based on completely 
different principles to directly search for the MLME estimator within the convex 
set in Section [3] This scheme couples two separate optimization techniques - 
scmidefinitc programming (SDP) and a derivative-free optimization method — to 
maximize the entropy over linear combinations of a maximal set of linearly independent 
ML estimators that spans the ML convex set. It will be shown that, owing to the 
mechanisms of SDP, one can make use of this scheme to infer the expectation values of 
a set of entanglement witnesses to verify the presence of entanglement in the unknown 
quantum state for composite systems. Finally, in Section 21 we establish a more 
robust numerical scheme that systematically generates the maximal set of linearly 
independent ML estimators that defines the convex set without fail. Instead of SDP, 
this scheme utilizes a nonlinear optimization routine that finds the global maximum 
of a highly nonlinear functional that is used to generate this maximal set. 

2. Brief review 

The likelihood functional \ogC({nj}; p) for a set of measurement data {nj} collected 
with a POM . IT, = 1 is given by 



where rij refers to the number of occurrences of the outcome IT, and pj = trjpllj}. 
The corresponding frequencies are given by fj = rij/N. In ML state estimation, 
the concave log-likelihood functional log C({nj}; p) is maximized to obtain the ML 
estimator /5ml 111 for the given set of data. If the number of linearly independent 
outcomes is D 2 , with D being the dimension of the Hilbcrt space, the data is 
informationally complete and the estimator /5ml is unique. In the case where this 
number is less than D 2 , the data is informationally incomplete and there exists now 
a continuous set of /5ml that yield the same maximal likelihood. This set is convex 
since the likelihood functional is concave in p. 

To choose the estimator /5mlme that has the highest entropy out of this convex 
set, we can consider the Lagrange functional 



(i) 



j 



X(\;p) = X(S(p) - S max ) + _l g£({n j };p) ) 



(2) 



X The symbol " " " is used to denote all estimators. 
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where A is the Lagrange multiplier corresponding to the constraint of maximal entropy 
S(Pmlme) = — tr{pMLME log Pmlme} = S max . In doing so, we maximize the two 
functionals S(p) and log(C({nj}; p)) simultaneously. We define pi a to be the estimator 
that maximizes I(A; p). 

If A = 0, I(X; p) = log C({nj}; p)/N and maximizing this Lagrange functional is 
just the procedure of ML. Since the data is informationally incomplete, there exists 
a convex plateau structure for the log-likelihood functional and maximizing I(A; p) 
yields a convex set of estimators. For large A values, the term XS(p) dominates, so 
that the resulting estimator pi t \-^-oo = 1/D. When A takes on a very small positive 
value [7J, the contribution from XS(p) becomes relatively much smaller than that of 
\og(C({nj}; p))/N, and any variation of the von Neumann entropy functional is only 
detectable over the state space region that coincides with the likelihood plateau. In 
other words, maximizing 1(0 < A — > 0; p) is equivalent to maximizing the entropy over 
the plateau. Therefore, />i,o<a->o = Pmlme- The iterative algorithm that is based on 
the steepest-ascent method is described in [7j-[9]. 

In practice, there is a limit to how small A can be. In particular, when A is 
smaller than some numerical threshold A t hres > 0, the gradient of XS(p) is no longer 
visible. In this case, the algorithm treats I(A < Athres; p) as I(A = 0; p) and performs 
ML estimation. Hence, the optimal parameter A is to be slightly above A t hrcs so 
that Pi t \>\ thrcs ~ pmlme- This inevitably introduces a small bias to the estimator. 
Determining the value of Athrcs analytically is quite complicated because Athrcs is 
actually a function of the true state and the POM, that is Athres = Athres({Hj }; ptrue)- 
Since ptmc is unknown to us, one needs to estimate Athres through repeated runs of 
the algorithm. In the next section, we will introduce an alternative scheme to search 
for the exact MLME estimator. 

3. Algorithm for entanglement and state verification 

To obtain the unique MLME estimator that has the highest entropy, a search has to 
be performed within the ML convex set. In this section, we introduce a numerical 
procedure to estimate pmlme in two steps. The first step is to obtain a collection of 
boundary states of the convex set for the measurement data obtained. In the next step, 
the operator /5mlme can be estimated using these boundary states with a standard 
function optimization routine. 

To carry out the first step, we need to identify the boundary of the ML convex 
set. Especially for large dimensions, the boundary of the convex set has an extremely 
complicated geometry that is too difficult to be analytically determined. Instead, we 
investigate its boundary by numerical means. We begin with the fact that the real 
functional f(H;p) = tr{pH}, where the operator H = W, is a linear functional of 
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p. Therefore, if we try to maximize (or minimize) f(H; p) over some subspace of p 
that has a well defined boundary, the maxima (or minima) of this linear functional 
are always on the boundary of this subspace. We can thus generate boundary states 
by maximizing or minimizing f(H; p), for a given Hcrmitian operator H, over the ML 
convex set. This problem is equivalent to the following optimization task: 

Maximize or minimize f(H;p) = ti{pH}, subjected to the 
following constraints: 

• p> 0, 

• tr{p} = 1, 

Ij} — L 



• tr{pILj} = tr{pMLlT, } for all j. 



This is a standard linear optimization problem, with linear and positivity constraints, 
that can be solved with the help of SDP [TU] . 

At this point, we need to find out the minimum number of boundary states that 
is required to search for /Smlme- To do this, we represent a generic ML estimator pml 
by a set of D 2 linearly independent tracc-orthonormal Hcrmitian basis operators 
Tj = rt satisfying the condition tr{r\,Tfc} = Sjk- With these basis operators, we can 
separate 

Anew Ainmeaa 

/3ml = £ a J i 7 eas + E ( 3 ) 

3=1 J=l 



— Pmeas Ptmmeas 

into the part in the measurement subspace of dimension D meas , and the rest of the 
state space that constitutes the unmeasured parameters (the unmeasured subspace) 
of dimension Ammeas, where D meas + D unmeas = D 2 , a 3 = tr{p M L r™ ^} and 
bj = trj^ML r" nmeas }. The generation of the basis operators r™ cas that span 
the measurement subspace can be done by the Gram-Schmidt orthonormalization 
procedure on the IIj-s, with all conventional inner products for vectors replaced by 
trace inner products for operators. In other words, the number of linearly independent 
POM outcomes lb, is D mea s- To generate the rest of the basis operators that span 
the unmeasured subspace, we continue the Gram-Schmidt procedure using randomly 
generated positive operators instead of the ITj-s. 

From ([3]) , it can be deduced that the maximum dimension of the ML convex set is 
-Dunmeas- To show this, we note that every ML estimator contains the same operator 
pmeas since the probabilities of pj = tr{p mcas Ilj} are fixed and trlpunmeaglL,} = 
for all IIj- by definition of the trace-orthonormal basis operators. This implies that 

§ If L Hcrmitian operators Aj arc all linearly independent, the rank of the matrix M with elements 
M jk = tr{A 3 A k } is L. 
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the only difference between any two ML estimators in the convex set is /5 unmeas . As 
the operators /5 moa s and /3 unmea s are linearly independent, it follows that any ML 
estimator can always be expressed as a linear combination of the unique p me as and 
the Dunmeas linearly independent basis operators that define p unmeas . This means that 
a set of .Dunmcas + 1 linearly independent boundary states is enough to look for the 
MLME estimator. As the operator /5 moa s is fixed by the measurement operators, the 
maximal number of free parameters that span the convex set is -Dunmcas- In some cases, 
however, the dimension of the ML convex set is lower due to the positivity constraint 
imposed on the ML estimators. In the extreme case, the convex set is restricted to a 
single point in state space. In these situations, we do not know its actual dimension 
and repeated generation of boundary states is necessary to estimate the maximum 
number of linearly independent boundary states. We remind the reader that with 
enough linearly independent states, the exact MLME estimator can be obtained up to 
numerical precision. In Section^ a different numerical procedure will be introduced 
to generate the maximal set of linearly independent ML estimators that spans the ML 
convex set without requiring any knowledge of the convex set. 

Apart from serving as a routine to numerically compute the boundary states of 
the convex set, SDP provides an additional useful function. For composite quantum 
systems, if one selects the Hcrmitian operators H to be entanglement witnesses W, 
one can obtain information about the presence of entanglement in the unknown true 
state even with informationally incomplete data. These witnesses have the properties 
that ti{p sop W} > for all separable states p sep and ti{p ent W} < for at least one 
entangled state p on t- The SDP routine, described above, thus looks for the maximum 
(or minimum) value of f(W; p) for any chosen witness operator W over the space of 
positive ps. This way, we can in fact infer a set of maximum (or minimum) witness 
expectation values over the ML convex set from this optimization procedure. The set 
of inferred maximum witness expectation values is particularly informative, for if the 
maximum value of f(W;p) for at least one of the randomly generated operators W 
is negative, we can immediately conclude that the true state is entangled since this 
state must lie within the ML convex set that results from the incomplete data. For 
practical computation, we can choose the witness operators W to be of the decom- 
posable form W = Q Tj , where Q is a positive operator with no product kets in its 
range and the symbol "t/' denotes a partial transposition with respect to the jth 
subsystem. For bipartite systems, these operators are optimal witnesses [HIE], that 
is no other witnesses can detect all entangled states that are detected by this witness, 
as well as other entangled states. For multipartite systems, these operators still serve 
as entanglement witnesses since tr{p sep Q T2 } = trjp^Q} > 0, although they are 
no longer optimal in general. To obtain a random set of boundary states of the ML 
convex set, random statistical operators Q are generated using the relation 
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Q = ^X}> < 4 > 
where X is a random operator which, when expressed in the computational basis, 
has complex matrix elements that are distributed according to the standard normal 
distribution of zero mean and unit variance. 

The second step involves an optimization procedure to maximize the entropy S(p) 

using the generated set of Mq < Z? U nmoas + 1 linearly independent boundary states 
f (i) ~l M ° 

| Pbd J ■ F° r the purpose of entanglement detection, these boundary states are ob- 
tained by maximizing the linear functionals f(W; p) of a set of witness operators W 
over the ML convex set according to the recipe described above. We start by writing 
a generic ML estimator as a linear combination of p^ inasmuch as 

M 

PML({^})=E^bd' (5) 
3=1 

where the tjS are normalized coefficients, such that J^j tj = 1, that are in general real 
such that Pml ({tj}) > 0. The task now is to look for the values of tj = £j nax for which 
the function S ({i" lax }) = S (pml ({^j nax })) is maximum over all real normalized coef- 
ficients. The unconstrained optimization of this Mo-dimensional function with respect 
to tj can be performed with any efficient multi-dimensional unconstrained optimiza- 
tion routine that is included in the standard libraries of commercialized mathematical 
softwares. In MATLAB, for instance, the function fminsearch does the job using 
the Nelder-Mead simplex method (NMS). When Mq is large, it is suggested in [13] 
that an adaptive version of the Nelder-Mead simplex method (ANMS) may be more 
advantageous in terms of shorter computation time. We take the resulting operator 
Psdp = Pml ({^j nax }) as the SDP MLME estimator. There is, however, a caveat to 
this optimization. Since the positivity of pml ({tj} ) is no longer guaranteed over the 
entire space of real normalized vectors, the entropy 



S({v 3 }) = -^VjlogVj , (6) 

3=1 

expressed in terms of the eigenvalues Vj of the statistical operator pml ({tj}) = 
\ v j) v j ( v j\ that is represented by its eigenbasis {|^)}, can take complex values. In 
order to restrict the optimization to yield only positive SDP MLME estimators, one 
can replace the entropy function in (j6|) with the conditional function 

f D 

-/jVj log Vj for puh ({tj}) > , 

S C ond({Vj}) = { j=l (7) 



k Sq < otherwise. 
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which effectively restricts the original search region to the admissible state space. 
To check if the evaluated operator /3ml ({tj}) is positive, a highly efficient way is to 
determine whether or not it admits a Cholesky decomposition. 

For full-rank SDP MLME estimators Psdp, there is a simple way to check that 
Psbp is indeed the MLME estimator. We recall that the form of the MLME estimator 
is given by 

PMLME ~ tr{e£*** n *} ' { ) 

where each Xj is a Lagrange multiplier for the constraint trjpMLME IL, } = trjpML IL,} 
for any ML estimator pml in the convex set. If /Smlme is full-rank, it follows from 
(O that the operator log /5mlme is a linear combination of only the POM outcomes, 
that is log pmlme resides in the measurement subspace. This is equivalent to the 
set of conditions trjFj 11 ™ ^ log/5 SDP } = for all r^ mea3 s. Defining the variables 
cj = tr{L" nmcas log /5sdp}, the quantity 7 = \J^2j c"j can be used to determine if 
Psdp is close to the actual MLME estimator up to some numerical precision. 

In summary, both the MLME estimator and information about the entanglement 
of ptruo can be obtained with informationally incomplete data using the following al- 
gorithm, which we coin as SDP MLME: 



SDP MLME 

(i) ML — Obtain the ML probabilities for the POM used to 
collect the measurement data with the ML algorithm. 

(ii) SDP — Perform SDP, using ((4]), to obtain the maximal set of 
Mo < -Dunmcas + 1 linearly independent boundary states that 
define the ML convex set using witness operators. Inspect the 
list of maximum expectation values of the witness operators 
and conclude that the unknown quantum state is entangled 
if any one of them is negative. 

(iii) Entropy maximization — Carry out function maximiza- 
tion with NMS or ANMS on the entropy function S ({tj}) 
with respect to all real normalized coefficients tj using ([5]) 
and to obtain /5sdp- 



Figured] summarizes the results. In Figure ([Taj) , the convergence of 7 is consistent with 
the fact that, in principle, the maximal number of nine linearly independent boundary 
states is enough to express the MLME estimator. The entanglement detection ratio in 
Figure [lb] is computed as the ratio of the number of detected pure states to the total 
number of generated states for every number of boundary states used. Ideally, this 
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Figure 1: Three plots showing various aspects of SDP MLME conducted on 100 random two- 
qubit entangled pure states as functions of the number of linearly independent ML boundary 
states used, with plots (a) and (c) averaged over all these pure states. A randomly-generated 
eight-outcome informationally incomplete POM is used throughout the simulations. In 
this case, every corresponding ML convex set is specified by nine linearly independent ML 
estimators. Decomposable optimal witness operators of the form Q^ 2 , where the operators 
Q are entangled pure states, are used to generate the boundary states with SDP. 
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ratio eventually goes to one if enough witness operators are generated to detect all the 
random pure states since there are no positive partial-transpose entangled states in 
this case [TTJ[T2]. For benchmarking, Figure [Tc] is generated to confirm the consistency 
of SDP MLME with the MLME algorithm described in Section [2] There exists an 
average bias in [Tc] that arises from a fixed Athres = 10~ 5 for all the pure states. 

4. Robust algorithm for incomplete state estimation 

Despite the usefulness of SDP in entanglement verification and boundary states 
generation as discussed in Section[3J the speed of the SDP routine strongly depends on 
the dimension of the Hilbcrt space and the total number of linear constraints imposed 
by the measurement data. When the total number of POM outcomes is large, there 
will generally be a considerable slowdown of the SDP routine as the search accounts 
for a large set of linear constraints in addition to the positivity constraint. Another 
feature of the SDP routine is that the sequence of boundary states that are generated 
from random Hcrmitian operators are not guaranteed to be linearly independent of 
one another. This means that typically, one would need to generate a large set of 
ML estimators that contains the maximal number of linearly independent estimators 
that span the ML convex set. For convex sets of large dimensions, this approach can 
be time-consuming. In this section, we propose a different search routine, in place of 
SDP, to directly look for linearly independent ML estimators within the ML convex 
set in a deterministic way. With this routine, we establish a feasible algorithm to look 
for the exact MLME estimator. 

To begin, we recall that a given set of Mq linearly independent ML estimators, 
as in implies the existence of a full-rank Mq x Mq positive Gram matrix Mq with 
elements given by 



This hints a straightforward strategy to cumulatively obtain the maximal set of linearly 
independent ML estimators that spans the entire ML convex set: Starting with a 
single ML estimator, the next estimator containing the same /5 mea s should be chosen 
such that the smallest eigenvalue cr m i n (Mq) of the Gram matrix Mq for these two 
estimators is maximized, and so forth, with the maximization performed over positive 

(i) 

estimators > 0. 

In general, a m i n (MQ) is a nonlinear functional of Mq that has multiple local sta- 
tionary points. This functional is also not differcntiablc and has undefined gradients 
at the boundary of the state space. To search for its global maximum, an appropriate 




(9) 
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numerical method to use is a nonlinear optimization algorithm that invokes pattern 
searches [14] and can cope with functionals that have ill-defined gradients. The solver 
for this algorithm, patternsearch, is readily available in MATLAB. The positivity 
constraint p^_. > is incorporated into the optimization algorithm using the aug- 
mented Lagrangian method with Cholesky decomposition. Another versatile feature 
of the proposed routine is that it can be set to terminate when the maximal number 
of linearly independent ML estimators is generated, such that the next ML estimator 
always yields a zero eigenvalue for Mq. We can, therefore, deterministically obtain 
the maximal set of linearly independent ML estimators that spans the ML convex set 
in this manner, in contrast with the SDP routine in SDP MLME. In this sense, the 
routine is operationally robust even without any prior information about the convex 
set. To ensure that the search is numerically stable, it is favorable to start the pattern 
search algorithm from a highly-mixed ML estimator. This is obtained by performing 
SDP as a couple of times and defining the starting estimator as an equal mixture of 
the resulting ML boundary estimators and the fairly mixed ML estimator obtained 
from the ML algorithm starting from the maximally-mixed state in computing the ML 
probabilities. Using this maximal set, the MLME estimator can be directly obtained 
by maximizing the entropy function in Q over all linear combinations of the linearly 
independent ML estimators in the set. These lead to the following pattern search 
MLME algorithm (PS MLME) for incomplete quantum state estimation: 



Verification of state and entanglement with incomplete tomography 



12 



PS MLME 

(i) ML ■ Obtain the ML probabilities, as well as the 
corresponding ML estimator that is fairly mixed, for the POM 
used to collect the measurement data with the ML algorithm 
starting from the maximally-mixed state. 

(ii) Carry out SDP a couple of times to obtain a few ML 
boundary estimators and check if the ML estimators are close 
to each other. If the average pairwise norm is smaller than a 
certain threshold, this means that the ML convex set can be 
approximated to be a single point and the ML estimator is 
taken to be the unique estimator. Otherwise, proceed to the 
next step. 

(hi) Definition of the ML convex set — Generate the maximal 
set of linearly independent ML estimators that defines the 
ML convex set by maximizing cr m i n (-MG) via the augmented 
Lagrangian pattern search algorithm. 

(iv) Entropy maximization — Carry out function maximiza- 
tion with NMS or ANMS on the entropy function S({tj}) 
with respect to all real normalized coefficients tj using ([5]) 
and (O to obtain the MLME estimator. 

The performance of PS MLME depends not only on the dimension L> U nmcas of the 
unmeasured subspace, but also on the complexity of the functional <7 m i n (Mq). More 
generally, as the dimension of the Hilbcrt space, or that of the unmeasured subspace, 
increases, the number of local maxima of er m m (Mq) increases. This translates to 
a longer computation time to locate the global maximum in the search for linearly 
independent ML estimators. Thus, for very large dimensions, PS MLME becomes 
inefficient and the approximate MLME algorithm, which is based on steepest ascent [7] , 
is a more practical substitute. On the other hand, PS MLME consistently gives more 
accurate MLME estimators as compared to the steepest-ascent algorithm, which yields 
biased results for large dimensions because of the finite A parameter. Hence, there 
is a tradeoff between computation time and the accuracy of the MLME estimators. 
Figure [2] compares the performances of the PS MLME and the steepest-ascent MLME 
algorithms for varying Hilbert space dimensions. Figure [2b] shows the consistently 
more accurate results obtained with PS MLME as compared to the steepest-ascent 
algorithm, in exchange for its longer computation time, illustrated in Figure I2al for 
higher dimensions. The simulations show that for the moderately large dimensions 
considered in Figure [2j much more accurate MLME estimators can be obtained using 
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Figure 2: Plots showing the average computation time for PS MLME (A) and steepest- 
ascent MLME (□) with Attires = 10 , as well as their respective average accuracies of the 
estimators for various Hilbert space dimensions D. The random POM for each dimension 
contains [D 2 /2J outcomes that are linearly independent. Each data point is averaged over 
50 randomly generated full-rank true states of every dimension. The computation time is for 
MATLAB running on an Intel i7 2.67 GHz Quad Core personal computer, where all function 
tolerances are set to 10~ 8 . 



PS MLME with longer computation time that is of the same order as that with 
steepest ascent. To reduce the computation time of PS MLME, one can consider an 
approximate maximization of <r m i n (Mq) as long as the resulting value is sufficiently 
large. Throughout the simulations, the duration of searching for each optimal ML 
estimator is restricted to five seconds. For even larger dimensions, PS MLME can be 
computationally demanding even when approximate maximization is carried out, and 
the steepest-ascent algorithm turns out to be a more realistic option. Efforts to further 
improve the performance of the PS MLME scheme are in the works. Nevertheless, we 
hope that the current work can serve as a stepping stone that helps to spur interesting 
discussions and novel contributions in related fields on numerical optimization over 
convex sets of positive operators. 
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5. Conclusion 

We have introduced a scheme to look for the unique estimator that maximizes 
the likelihood and entropy for informationally incomplete measurement data. This 
involves two main procedures: a generation of linearly independent boundary 
maximum-likelihood estimators that spans the convex set with semidefinitc 
programming, and an entropy maximization with these estimators using a standard 
function optimization routine. Furthermore, for composite quantum systems, one can 
apply this scheme to infer the expectation values of a set of entanglement witnesses. 
This information allows us to verify the entanglement of any quantum state with 
informationally incomplete data. However, semidefinitc programming does not offer a 
definite control over the generation of maximum-likelihood estimators. This motivated 
us to develop an alternative scheme that is more operationally robust than the former 
one to search for the maximum- likelihood- maximum-entropy estimator. This latter 
scheme makes use of the pattern search optimization algorithm that is suitable for 
maximizing a nonlinear function required to deterministically generate the maximal set 
of estimators that defines the convex set. With numerical simulations, we showed that 
the latter scheme gives much more reliable results than the MLME algorithm discussed 
in Section [5] at the expense of slightly longer computation time when the dimension 
of the reconstruction Hilbert space, or the unmeasured subspace, is moderately large. 
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